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We study numerically the fully nonlinear gravitational collapse of a self-gravitating, minimally- 
coupled, massless scalar field in spherical symmetry. Our numerical code is based on double-null 
coordinates and on free evolution of the metric functions: The evolution equations are integrated 
numerically, whereas the constraint equations are only monitored. The numerical code is stable 
(unlike recent claims) and second-order accurate. We use this code to study the late-time asymp- 
totic behavior at fixed r (outside the black hole), along the event horizon, and along future null 
infinity. In all three asymptotic regions we find that, after the decay of the quasi-normal modes, 
the perturbations are dominated by inverse power-law tails. The corresponding power indices agree 
with the integer values predicted by linearized theory. We also study the case of a charged black 
hole nonlinearly perturbed by a (neutral) self-gravitating scalar field, and find the same type of 
behavior — i.e., quasi- normal modes followed by inverse power- law tails, with the same indices as in 
the uncharged case. 



PACS number(s): 04.70.Bw, 04. 25. Dm 
I. INTRODUCTION 

The no-hair theorems state that except for the mass, 
the electric charge, and the angular momentum, all the 
features of fields which collapse to a black hole will be 
unobservable to external observers at late times. It is 
therefore interesting to study the mechanism by which 
the hair is radiated away (or absorbed by the black hole). 

Until recently, the late-time evolution of non-spherical 
gravitational collapse was investigated primarily in the 
context of linear theory. That is, the deviations from 
spherical symmetry were considered as infinitesimally- 
small perturbations over a fixed curved background. The 
late-time behavior of such perturbations has been studied 
for three different asymptotic regions: (a) at fixed r, (b) 
along null infinity, and (c) along the future event hori- 
zon (when the collapse is to a black hole) . Qualitatively, 
the evolution of the linearized perturbations is similar in 
these three asymptotic regions: During the first stage, 
the perturbations' shape depends strongly on the shape 
of the initial data. This stage is followed by the stage of 
quasi-normal (QN) ringing, in which the perturbations 
oscillate with an exponentially-decaying amplitude. The 
corresponding complex frequency is characteristic of the 
parameters of the background black-hole, and is indepen- 
dent of the details of the initial perturbation. Finally, 
there are also 'tails', characterized by an inverse power- 
law decay. 

The asymptotic region (a) was first studied by Price 
|jlj, who analyzed the linear perturbations over a fixed 
Schwarzschild background. Price found that after the QN 
ringings die out, the perturbations at fixed r (outside the 
black hole) decay according to t~( 2l+fJ,+i \ w here fx = 1 
if there were an initial static mode, and /i = 2 otherwise. 
Here, I is the multipole moment of the mode in question, 



and t is the standard external Schwarzschild time coor- 
dinate. Asymptotic regions (b) and (c) were considered 
by Gundlach, Price and Pullin p|, who showed that the 
tails along null infinity decay according to u e where 
u e is the outgoing Eddington-Finkelstein null coordinate. 
(Hereafter we use the notation u e and v e for the outgo- 
ing and ingoing Eddington-Finkelstein null coordinates, 
correspondingly, in order to distinguish them from other 
types of null coordinates which we use later.) Along the 
event horizon (EH), the inverse-power indices were found 
to be similar to the asymptotic limit (a); Namely, the 
'tails' decay according to w c T^ 2i+p+1 ' > ■ 

Recently, Krivan, Laguna and Papadopoulos studied 
numerically the evolution of linearized scalar-field || and 
spin-2 perturbations Q over a fixed Kerr background. 
(See also ||.) They concluded that 'tails' are expected 
also for the Kerr background, with power-law indices sim- 
ilar to those obtained for Schwarzschild. This provides 
additional motivation for the study of the fully-nonlinear 
evolution of perturbations in the spherically-symmetric 
case as a toy-model for the spinning case, since the spher- 
ical case is much simpler to deal with (both analytically 
and numerically). 

The numerical simulation of the fully-nonlinear 
gravitational collapse of a spherically-symmetric self- 
gravitating scalar field was recently carried out by 
two groups: Gundlach, Price and Pullin || (GPP), 
and Marsa and Choptuik (MC) @. In both anal- 
yses, the coordinates used were non-vacuum general- 
izations of the (one null+r) outgoing |6) or ingoing 
[Q Eddington-Finkelstein coordinates. These numerical 
analyses demonstrated the QN ringing as well as the 
power-law 'tails' for lines r = const. In addition, MC 
also demonstrated the power-law decay at the EH. 

In this paper, too, we study the nonlinear spherical 
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gravitational collapse of a self-gravitating scalar field. 
However, we shall use different coordinates, different nu- 
merical methods, and a somewhat different model. Our 
numerical code is stable and second-order accurate and is 
based on free evolution and double-null coordinates. This 
combination has several advantages: First, the null coor- 
dinates are very well adapted to the hyperbolic character 
of the field equations evolved: The evolution is along the 
characteristics, and consequently there is no restriction 
analogous to the Courant condition. Second, in double- 
null coordinates the interpretation of the causal structure 
of the numerically-produced spacetime is trivial. Also, 
double-null coordinates can be chosen such that the met- 
ric is regular at the EH, which is not the case in outgoing 
Eddington-Finkelstein coordinates . 

Our coordinates and integration scheme allow us to 
study the evolution to arbitrarily late times, and there 
is no need to introduce an artificial outer boundary |Q. 
Our analysis demonstrates both the QN ringing and the 
power-law tails. One of our main objectives in this in- 
vestigation is to numerically determine the power-law in- 
dices of the late-time tails in the nonlinear collapse prob- 
lem, and to compare them to the predictions of the lin- 
ear perturbation theory. In cases (a) and (c), we ob- 
tained power-law indices similar to those found by ||f7| , 
though with improved accuracy. In addition, we ob- 
tain the power-law index at null infinity, which was not 
studied so far in nonlinear collapse. In all three asymp- 
totic limits, we find an excellent agreement between our 
numerically-obtained indices and the values predicted by 
the linear perturbation analyses 0J^]. (Such an agree- 
ment is expected, even in a very nonlinear collapse prob- 
lem, because of the "no-hair" principle — see e.g. Ref. 

Whereas this paper considers only the external part of 
the black hole, we are currently investigating the inner 
structure of charged black holes with a similar numerical 
code B. In fact, our main motivation in this project is 
to develop the numerical approach and techniques which 
could later be used in investigating the black hole's inte- 
rior. The determination of the correct late-time power- 
law index is essential for that purpose. 

Similar self-gravitating collapse scenarios have been re- 
cently used for the study of critical phenomena in black 
hole formation [p|~[i"T||. At its present form our code is 
incapable of treating these phenomena because we have 
not attempted to include the neighborhood of the ori- 
gin in the domain of integration. The configurations we 
are interested in here are by far super-critical, and the 
aspects which concern us do not require the integration 
near the origin (see below). 

This paper is organized as follows: In Section II we 
present the model for the collapse, and the correspond- 
ing field equations. Section III describes the numerical 
approach, and Section IV discusses the stability and ac- 
curacy of our code, and the tests used to verify them. It 



has been recently argued |T^] that unconstrained codes in 
double-null coordinates suffer from inherent instabilities. 
We show that this is not the case, and that our code is 
indeed stable and converges with second order. (In the 
Appendix we explain this in greater detail.) In Section 

V we present our numerical results for the collapse of a 
scalar field over a Minkowski background, leading to the 
formation of a Schwarzschild black hole, and in Section 

VI we consider the collapse of a (self-gravitating, neutral) 
scalar field on a RN background. Finally, in section VII 
we summarize and discuss our results. 



II. THE COLLAPSE MODEL 
A. The field equations 

We shall consider the spherically-symmetric gravita- 
tional collapse of a self-gravitating, minimally-coupled, 
massless scalar field. In the uncharged case, the system 
is described by the coupled Einstein-Klein-Gordon field 
equations. We shall also consider the charged case, i.e. 
the case in which a (source-less) spherically-symmetric 
electric field is also present. In this case, the system is 
described by the coupled Einstein-Maxwell-Klein-Gordon 
field equations. 

We write the field equations in double-null coordinates. 
The line element takes the form 



ds 2 = —f(u, v) du dv + r 2 (u, v) dVL 2 



(1) 



where dfl 2 is the line element on the unit two- 
sphere. The general spherically-symmetric solution of 
the Maxwell equations in these coordinates is 
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otherwise, where Q is a free parameter, 
which is interpreted as the electric charge, and where 
Ffj, v is the Maxwell field tensor. The contribution of this 
electric field to the energy-momentum tensor is 
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The energy-momentum tensor of a massless scalar field 
$ is 



(4) 



This field satisfies the Klein-Gordon equation <F Q ; " = 0, 
which, in our coordinates, takes the form 
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uv-r- {r,u®,v +r, v $, u ) = 0. 
r 



(5) 
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The Einstein field equations are = 87rT M „, where 
the energy-momentum tensor is the sum of the contri- 
butions of both the electromagnetic and scalar fields, 
These equations include two evolu- 
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tion equations 



and 
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f,uf,V , j. f 1 r. 



(6) 
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supplemented by two constraint equations: 
r,uu - (ln/) )U r,„ + r($ u ) 2 = 

?>„ - (ln/) )V r,„ + r($ )V ) 2 = 0. 



(7) 



(8) 



The constraint equations are not independent of the dy- 
namical equations: Any solution of the evolution equa- 
tions will also be a solution of the constraint equations, 
provided only that the latter are satisfied on the initial 
hypersurface. (This is assured by virtue of the contracted 
Bianchi indcntities.) 



B. The formulation of the characteristic problem 

In our numerical scheme, we shall use the three 
equations @— (@) to evolve the three unknowns r(u,v), 
f{u,v), and $(«,«). These equations form a hyperbolic 
system, and thus ensure a well-posed initial-value for- 
mulation. In the double-null coordinates we use, it is 
most natural to use the characteristic initial-value for- 
mulation, in which the initial values of the unknowns 
(but not of their derivatives) are specified on two null 
segments, u = const = m and v = const = Uj. 

In such a numerical scheme (often called free evolu- 
tion), the constraint equations are only imposed on the 
initial hypersurfaces. As mentioned above, the consis- 
tency of the evolving fields with the constraint equations 
is mathematically guaranteed. We use the constraint 
equations to check the accuracy of the numerical sim- 
ulation. 

From the pure initial- value viewpoint, we need to spec- 
ify three initial functions on each segment of the initial 
surface: r, /, and $. The constraint equations, how- 
ever, reduce this number by one: Eq. (|) imposes one 
constraint on the initial data at v = i>i, and similarly, 
Eq. (||) imposes one constraint on the initial data at 
u = m. The remaining two initial functions, however, 
represent only one physical degree of freedom: The other 



degree of freedom expresses nothing but the gauge free- 
dom associated with the arbitrary coordinate transfor- 
mation u — > u(u) , v — > v(v) (the line element ([j]) and all 
the above equations are invariant to this transformation) . 
In what follows we shall use a standard gauge, in which 
r is linear with v or u, correspondingly, on the two initial 
null segments. On the outgoing segment we take r„ = 1. 
On the ingoing segment, we take r u — const = r u o- The 
initial values of r are thus uniquely determined by the 
parameter tq = r(ui,Vi). We choose Ui = and Vi = ro, 
and thus we find: 



r u (u) = r + ur u0 



(Hereafter, we denote the initial values of the three fields 
on the two initial segments by r u (u), f u (u),<f> u (u) and 
r v{v), f v (v), $v(v), correspondingly.) Then, we can freely 
specify $«(«) and $> v (v) (this choice represents a true 
physical degree of freedom) . The initial value of / is now 
determined from the constraint equations, namely 



(9) (In f u ),u = r„ ($„,«) 2 /r„o , (ln/«) ($„,,) 2 , (10) 



together with the choice f(ui, Vi) — 1. Thus, in the gauge 
we use, the geometry in the entire domain of dependence 
is uniquely determined by the two initial functions $«(«) 
and $ v (v), and the two parameters ro and r U Q. (Later we 
shall relate r U Q to the initial black- hole mass.) In what 
follows, we shall consider initial data corresponding to 
a compact ingoing scalar-field pulse, over a background 
of either Minkowski, Schwarzschild, or RN. Namely, we 
shall assume that $ u (u) = 0; and is also zero, except 
at a finite interval v± < v < V2 (with some v\ > v,-) where 
<I>t, ^ 0. For concreteness, in the range v\ < v < V2 we 
shall take = Asm 2 [TT(v —v\)/{v2 — v\)\ . This choice 
for the initial data is smooth at the matching points v — 
v\ and v = The determination of f v (v) throughout 
u = Ui, by analytically integrating the second constraint 
equation in (|l^), is straightforward for such a pulse. On 
the ingoing segmant v — Vi we have f u (u) = 1- 

The geometry is static (with $ = 0) in the entire range 
v < Vi, with a mass parameter Mq. To relate Mq to 
the above initial-value parameters, we define the mass- 
function M(u, v) HU by r )A4 r" = l-2M(u,v)/r+Q 2 /r 2 . 
In our coordinates this becomes 



M(u, v) = (r/2)(l + 4f-\ u r,v) + Q 2 /2r 



(11) 



which yields M = (r /2)(l + 4r„ ) + Q 2 /2r . Thus, 
our initial- value set-up is determined by the initial mass- 
parameter Mq, the charge Q, and the perturbation am- 
plitude A (together with the auxiliary parameters V\ , v-i 
and ro). 

We shall particularly study two cases: (i) M = 0, Q = 
(Minkowski background), and (ii) M = 1, Q ^ (RN 
background). [Note that no loss of generality is caused 
by the choice Mq = 1 , because of the scale- invariance na- 
ture of the problem.] We shall not elaborate here on the 
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situation of self-gravitating scalar field collapsing over a 
Schwarzschild background, as its outcome is qualitatively 
similar to case (i) (and, the collapse over Minkowski 
brings out the nonlinear aspects in a sharper way). We 
shall use, however, the pure (A — 0) Schwarzschild case 
as a test-bed for our numerical code. 

In what follows, we shall use the symbols u and v to 
denote the outgoing and ingoing null coordinates in the 
specific gauge described above. Notice that v is closely 
related to w e at u > M, and u is Kruskal-like near the 
EH (namely, it regularizes the metric function / at the 
EH). 

The double-null line element suffers from a non- 
physical coordinate singularity at the origin (i.e., the 
timelike worldline r — 0, where the geometry is perfectly 
regular). This singularity may cause difficulties in the 
numerical study of case (i) above (i.e., Minkowski back- 
ground). In order to overcome this difficulty, we restrict 
the domain of numerical integration in this case such that 
it will not include the origin. That is, the characteristic 
initial segment v — Vi ends before it reaches r = 0. Since 
it is very essential that the domain of integration will 
include the EH, we must demand that the ingoing ray 
v = Vi will intersect the EH before it intersects r = 0. 
This is achieved if the amplitude parameter A is suffi- 
ciently large. 

III. THE NUMERICAL CODE 

Our numerical code is based on the standard proce- 
dure for second-order integration of (1+1) hyperbolic 
equations in double-null coordinates: Let du and dv be 
the finite increments in the u and v directions, respec- 
tively. Let us also denote schematically the three un- 
knowns r, /, $ as hi , i = 1,3] These unknowns satisfy a 
field equation of the form 

hi : uv — -Fi(hj 7 hj U: hj v *j , j — 1,3 . (1^) 

Assume now that we already know the values of hi at the 
three grid points p\ = (uq,vq), p 2 = (uq + du,Vo), and 
P3 = {uq. vq + dv), and we would like to evaluate them 
at j»4 = (uq + du, vq + dv). We then use the substitution 

h^^-hf+hf^+hf+F^dudv , (13) 

where, for any function g, g^ = g(pk), and ps is the 
intermediate point: p§ = (u + du/2,VQ + dv/2). In order 
to evaluate the functions hj, hf^, h°l (required for the 

determination of ) to the desired accuracy, we use the 
standard "predictor-corrector" method. This procedure 
results in a second-order accuracy. With this method, 
we first calculate hi along the ingoing ray v = Vi + dv — 
starting at u = Ui + du, then solving for u — Ui + 2 du, 
and so on, until the last grid point at u = Uf, Then we 



turn to the next ingoing ray v — Vi + 2 dv, and solve for 
all grid points along this ray; By this way we solve, ray 
after ray, until we cover the entire domain of integration, 
Vi < v < Vf , Ui < u < Uf. The accuracy is controlled 
by the global grid parameter N, which is the (initial) 
number of points per a unit interval in both the v and u 
directions. Typically we used N = 10, 20 or 40, though 
in certain cases we also used the values 5, 80, and 160. 

As long as our domain of numerical integration does 
not include the EH, this numerical scheme can be used 
with a fixed grid without any difficulties. When the EH 
is included, however, we face a fundamental difficulty. 
For simplicity, assume at this stage that <& = and 
the spacetime is Schwarzschild (though the same con- 
ceptual difficulty arises also in non-static spacetimes). 
Let us denote by Uh the value of our Kruskal-like co- 
ordinate u at the EH. Let uq be a grid value of u just 
near the EH, and let u\ be the next grid point in u, i.e. 
ui = uq + du. We define Sr(v) = r(u\,v) — r(uo,v). Of 
course, for the validity of the numerical integration it is 
necessary that Sr <C r — and we shall indeed select the 
grid parameter du sufficiently small so as to satisfy this 
requirement at the initial segment v = Vi. The prob- 
lem is that, Sr grows unboundedly and very rapidly with 
v. Thus, in terms of the Eddington-Finkelstein coordi- 
nate v e , along the horizon Sr oc exp(i; e /4M) (because 
at the horizon of Schwarzschild, r u oc exp(w e /4M), and 
8r(v) = r u (uo, v) du.) It is therefore obvious that a code 
based on a fixed du cannot be used here. One might at- 
tempt to use a numerical scheme in which du depends on 
u (but not on v) in such a way that it becomes extremely 
small at the EH. But this turns out to be impractical 
too, because of the extremely large exponential factor: 
Typically we need to integrate up to v e values of at least 
a few hundreds times M (otherwise we cannot study the 
power- law tails with a sufficient accuracy). This would 
demand a value of du as small as, say, 10 near the 
EH, which is obviously impractical, due to the roundoff 
error and other reasons. 

In order to overcome this difficulty, we must use a 
dynamical grid-refinement algorithm. A sophisticated 
dynamical refinement scheme was recently developed by 
Hamade and Stewart [|o) , in order to analyze the critical 
behavior at the origin. For our purposes, however, it is 
sufficient to use a simpler refinement scheme, which we 
call point splitting: In certain values of v, we check the 
variations in r (and, in fact, in all hi) between any two 
adjacent grid points. If the difference in r between such 
two points pi = (uq, vq) and P2 = (uq + du, Vq) is greater 
than some threshold value, we introduce an intermedi- 
ate grid point p' 2 = (uq + du/2,Vo), and calculate the 
interpolated values of all unknowns hi at that point. We 
can now use the above three-points integration scheme to 
calculate hi at p' A = (uq + du/2,VQ+ dv) from the values 
of these fields at pi,P2,P3 (and, later, to calculate hi at 
Pi = (uq + du, vq + du) according to the field values at 
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P2)P2,P4)- This numerical procedure functions very well 
in double-null coordinates, especially due to the follow- 
ing reasons: First, in the three-points integration scheme, 
{PiiP2,Pz) — * Pa, there is no reference to any grid points 
at u < uq or u > uq + du. Therefore, it does not matter 
whether the increments du are uniform or not. Second, 
in this scheme there is no restrictions on the ratio of du 
and dv. 

In practice, we proceed as follows. We register all grid 
values of u (at a given v) in a vector u(I), where I is an 
integer index. The values of the three unknowns are reg- 
istered in corresponding three vectors hi{I). We define 
three threshold parameters h\ for the three unknowns, 
and also a "band parameter" Vb (typically we take Vb 
to be of order M). At the end of each interval Vb in 
v (a "band"), we check the variation of all three fields 
along the vectors hi(I). If for a given I the relative dif- 
ference \[hi(I + 1) — hi(I)]/hi(I)\ is found to be greater 
than hi (for any i), then we add a new grid point at 
u = [u(I) + u(I + l)]/2. In such a case, we calculate 
the values of the fields hi at that new point by inter- 
polation (usually we perform a four-points interpolation, 
based on, e.g., u(I — 1), u(I), u(I + 1), u(I + 2)). We now 
update the vectors u(I) and hi(I), by assigning the value 
/ + 1 to the new grid point. (Before creating this new 
grid point, we arrange an empty "slot" for it, by shift- 
ing all grid index values I' > I by one.) The threshold 
values hi are taken to be proportional to l/N, in order 
to preserve the rule of N as a parameter that controls 
the global accuracy (that is, the number of grid points 
in the u axis should be proportional to N). The band 
parameter Vb is taken to be independent of N. 

Because our goal is to study the evolution in the en- 
tire black-hole exterior, the domain of integration must 
include the EH and thus extend into the black hole (i.e., 
Uf > Uh). Then, if Q = 0, the numerical integration 
will terminate at some finite v, beyond which the ingo- 
ing null lines v — const intersect the spacelike r = 
singularity (before u — uf). In order to overcome this 
difficulty, we simply chop the vector u{I) just beyond 
the apparent horizon (AH). Namely, at the end of each 
band, we first find I ah, the value of the index / where 
the AH is located. This is the value of / satisfying 
r, v (u = u(I - 1)) > and r >v (u = u(I)) < 0.0 We then 
chop the vector u(I) at, say, I — I ah + 1- This ensures 
that the domain of integration never gets close to the 
spacelike r = singularity — and yet it contains the entire 
external part of the domain Vi < v < Vf , u,; < u < Uf, 



up to (and including) the EH. 

With these procedures of point-splitting and chopping, 
our code can in principle run to arbitrarily large v val- 
ues. Due to point-splitting, however, the number of grid 
points in the vector u(I) grows linearly with v e , so the 
integration time grows like v e 2 . In order to significantly 
decrease this time, we introduce two additional types of 
numerical manipulations: 

point removal: The successive addition of points near 
the event horizon results in a coverage of the regions 
r 3> 2M, v e 3> AI with an approximately uniform density 
of order N points per unit u e (Recall that a point added 
at u just before the EH will, after an interval v e 3> M, ap- 
proach r S> 2M.) However, for an appropriate coverage 
of the variations in r, a much smaller density of about 
NM/r points per unit u e will be sufficient (note that 
at r > 2M, r is approximately linear in u e along lines 
v = const). Thus, in order to save integration time, at 
the end of each band we check along the vector u(I) and 
simply remove all unnecessary points, thereby shortening 
this vector. The criterion for necessity or otherwise of a 
point I is qualitatively similar to that of point-splitting: 
Again, we define the threshold values h'l for point re- 
moval (typically, h'l is slightly smaller than hi/ 2). If for 
all i we have \{hi(I - 1) - hi(I)]/hi(I)\ < h'l, then the 
point I is removed^]. 

gauge correction: As it turns out, for any Kruskal- 
like u (which is necessary for a regular coverage of the 
EH) and Eddington-like v, f grows exponentially with v e 
along the EH. In addition, along lines v — const ^> M, 
f grows rapidly (exponentially with u e ) in most of the 
interval m < u < Uh- The above numerical scheme 
handles very well this behavior of /. However, the sig- 
nificant variation of / with u implies that points can 
hardly be removed, which results in a long computa- 
tion time. In order to overcome this difficulty, we in- 
troduce (as an option) a gauge correction at the end of 
each "band." That is, we perform a coordinate trans- 
formation u — ► u new (u). The value of / is gauged ac- 
cordingly: f new = (duoid/ du new )f id (the variable r is 
unchanged). Our field equations are invariant to such a 
coordinate transformation. The function u new (u) is to 
be chosen so as to decrease the variation of f new with 
u. A convenient choice is to take u new (u, vq) = r(u,Vo) 
(where vq is the value of v at the ingoing ray where 
the gauge transformation is carried out), in which case 
fnew turns out to be approximately constant through- 
out the ingoing ray. (Another convenient choice is to 



* We also calculate uah by interpolating between the two 
points u(I — 1) and u(I). (In the Schwarzschild or RN cases 
uah = Uh , but in the general dynamic case uah > Uh ■ Recall 
that only the AH can be found locally. However, for large v 
the AH should coincide with the EH.) 



'In fact, the criterion we use for the variation in <3? (for both 
point-splitting and point removal) is somewhat more involved: 
It refers to the variations in both $ and <I> iU . This is essen- 
tial for the appropriate coverage of the maxima and minima 
regions in the QN ringing. 
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define u new {u) by the demand f n ew{u,Vo) = 1.) We re- 
call, however, that our goal is to numerically compute 
f original (as well as t and $) as a function of 

^oriqinal 

and u, and the gauge transformations are just a sub- 
sidiary manipulation. In order to accomplish this goal, 
we must keep record of two additional variables: (i) 
the vector u or i g i na i(I) , and (ii) the vector R g (I), where 
R g = ducurrentl du or i g inai is the cumulative gauge fac- 
tor. At each gauge transformation, the latter is updated 
according to R g ( new ) = {du new / du old )R g{old) . Conse- 
quently, in the original gauge the metric function / is 
given by f original = R g /current- Hereafter, whenever 
we mention u and /, we refer to u origina i and f „ g i„ a i, 
correspondingly. 

The combination of the above four types of numerical 
manipulations yields an accurate and efficient numerical 
code. It is important to recall that, whereas the point- 
splitting and chopping are necessary for an integration to 
large values of v, the point removal and gauge correction 
are optional, and are aimed to save integration timeM. 
With these two manipulations, the typical number of grid 
points in the vector u(I) grows only logarithmically with 
v e , instead of linearly. In practice, the integration times 
in long runs are reduced by a factor of 10 or so. Typically, 
in a running to large v, almost all points in u are "born" 
in point-splitting near the EH, and are later removed 
when they approach r ^> 2M. We emphasize that in the 
numerical scheme described here the increment in v is 
fixed, dv = l/N. 

IV. STABILITY, ACCURACY, AND ERROR 
ANALYSIS 

Stability 

Gundlach and Pullin (GP) recently argued that 
any free-evolution scheme will be inherently unstable, in 



' In a point-splitting, the variables u orig i na i(I) and R g (I) 
are interpolated at the point added, like the other variables 
hi(I). 

S If a gauge-correction is not used for any reason, then, as 
a consequence of successive point-splittings, the difference in 
u between two adjacent points near the horizon becomes as 
small as, say, 10~ 100 . Then, due to the roundoff error, it is 
not possible to calculate du(I) = u(I + 1) — u(J) directly at 
each step. One way to overcome this difficulty is to keep an 
independent vector du(I), and to update it at every point- 
splitting by dividing du(I) by two. (Recall that it is du that 
is involved in the finite-difference integration scheme, not u.) 
Another possibility is to shift u by a constant, e.g. at the end 
of each "band," so as to assign the AH the value u — — in 
this case u(I + 1) — u(I) can be calculated directly at each 
stage. 



the sense that small violations of the constraint equa- 
tions will grow exponentially with t along lines r = const 
- even if the evolution equations are exactly satisfied. 
We disagree with the theoretical analysis and interpre- 
tation made by GP, for reasons explained in the Ap- 
pendix. Also, our numerical tests did not indicate any 
such numerical instability. Certain entities exhibit an 
exponential growth, but these entities are not the ones 
that may be used as authentic error indicators; rather, 
the exponential growth we encountered is merely a re- 
flection of the passive exponential growth exhibited by 
various gauge-dependent entities (e.g., r >u , for Kruskal- 
like u) along lines r = const (or along the EH) in the 
Schwarzschild geometry. We discuss this issue exten- 
sively in the Appendix. 

In our stability tests, we numerically reconstructed the 
Schwarzschild spacetime (as well as RN and other space- 
times) up to t values of many thousands times M, and 
with values of the grid-parameter N ranging from 5 to 
160. In all these cases, we found a stable numerical evo- 
lution. 

Accuracy checks and error analysis 

We used several methods to test and monitor the accu- 
racy of our numerical code: (i) Comparing the results ob- 
tained with different values of the grid-parameter N; (ii) 
monitoring the discrepancy in the two constraint equa- 
tions (||,[|) (as explained above, our integration scheme 
does not involve these equations); (iii) numerically repro- 
ducing known exact solutions: The Minkowski solution, 
the vacuum Schwarzschild solution, the electro-vacuum 
RN solution, and the self-similar spherically symmetric 
scalar-field solution [|l5| for the Einstein-scalar field equa- 
tions. In all these cases our code reproduced the exact 
solutions very well (we have compared the metric coeffi- 
cients r, / and the mass function) . All these tests indi- 
cated that the code is stable, and the numerical errors 
decrease like N~ 2 , as expected from a second-order code 
(see examples below). 

We present here the results for the numerical repro- 
duction of the Schwarzschild solution. (Similar results 
were obtained for the other above-mentioned tests.) We 
start with initial data corresponding to M = 1. The 
drift of M from its initial value may then be used as an 
error indicator. Figure ^ displays M as a function of t 
along a line r = const, Fig. || shows the drift of M as 
a function of v along the EH, and Fig. || shows M as a 
function of u e along null infinity, for various A~ values. 
From these Figures it is apparent that the mass drift is 
linear with time, and decreases like N~ 2 , as expected. 
We also compared the mass function obtained from lo- 
cal differentiation [Eq. (pT[)] with the dynamical mass 
function which evolves according to the wave equation 
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r 3 o 9 / 2m Q 2 \ 
m, uv = 2j<5> 2 u <$>l-r ( 1 - — + ^-J (14) 

Both expressions for the mass function agree with each 
other (in the limit of large N ) . 




200 400 600 800 1000 1200 1400 1600 1800 2000 
t 



FIG. 1. The drift of the mass function along r — const. 
Shown here is the mass function as a function of t for 
N = 10, 20 and 40. We took here r = 8. 



1.002 




100 200 300 400 500 600 

FIG. 2. The drift of the mass function along the EH, as a 
function of v, for N = 10, 20 and 40. 

Another check we performed was to compare the values 
of / as a function of t along lines r = const of the numer- 
ically reproduced Schwarzschild spacetime and the exact 
analytical counterpart |l^] . Figure f| displays the results 
we obtained for r = 3M. (Similar behavior was found for 
other values of r.) It is convenient to compare the values 
of / in the outgoing Kruskal coordinate Uk and the in- 
going Eddington coordinate v e . For the Schwarzschild 
solution one finds that g UkVc = iMl e -r/(2M) e v e /(AM) _ 
Figure [| shows the ratio F between gu k v e for the ex- 
act Schwarzschild solution and gu h v e m the numerically 
reproduced spacetime, along r — 3M, vs. the ingoing 
coordinate v, for several values of the grid parameter N. 



This figure clearly indicates the second-order convergence 
of the code to the correct theoretical value. 
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FIG. 3. The drift of the mass function along null infinity, 
as a function of u e (calibrated such that u e = on u = 0), 
for TV = 10, 20 and 40. 
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FIG. 4. The ratio F of gu k v e of the exact Schwarzschild 
solution and gu k v c of the numerically reproduced spacetime. 
Shown are the values for N — 40, 80, and 160. The conver- 
gence indicates a second-order code. The deviation of the 
curves from straight lines results primarily from the linear 
drift of M. 

As we mentioned above, we also use the constraint 
equations to monitor the errors. Let us denote the en- 
tities in the left-hand side of Eqs. d||J) by C" u and C' v , 
respectively. Now, as they stand, C' u and C' v cannot be 
used as measures of the intrinsic local error in the repro- 
duced spacetime, because they are not gauge invariant. 
Instead, whenever r „ or r „ are non-vanishing, we may 
define 

C u = C u lr 2 u , C v = C' v /r 2 v . 

Since C u and C v are gauge-invariant, they provide an 
invariant measure of the local numerical error. An altcr- 
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native gauge-invariant indicator is C = / _1 (C^C^,) 1 ^ 2 . 
Note that the indicator C v cannot be used at the AH, 
where r t vanishes. Instead, one may use the indicator C 
there. Figure ^ displays C u and C v at constant r as func- 
tions of t (we took here r = 3), and C u and C at the EH, 
as functions of v. From Figure || one can see that these 
indicators are roughly constant with time. (The noise is a 
results of the second-order numerical differentiation nec- 
essary for the computation of the indicators.) In partic- 
ular, no exponential growth occurs. This demonstrates 
the stability of the code. We found a similar behavior 
also for the electro-vacuum RN spacetime. 




100 200 300 400 500 600 

time 



FIG. 5. The gauge- invariant error indicators \C U \ and |C„| 
along a line r — const as functions of t, and \C U \ and \C\ along 
the EH as functions of v e . Line a: \C U \ along r — const, line 
b: \C U \ along the EH, line c: \C V \ along r = const, and line d: 
|C| along the EH. The data are taken for r — 3 and TV = 40. 

Figure ^ shows the rate of convergence of various er- 
ror indicators as N increases. The spacetime simulated 
here is RN with Q/M — 0.8. (The other exact solutions 
we checked produced similar results.) Shown are the l p - 
norms, for several p values, of the two vectors made of the 
values of the indicators C u and C v , respectively, along a 
particular outgoing ray located before the EH. We used 
the following values of N: 5,10,20,40,80, and 160. The 
apparently straight lines in the logarithmic graphs indi- 
cate a second-order convergence^. (The break in the 
lines e and f for N — 160 seems to be a roundoff effect.) 
The other error indicators we used (e.g., the drift of the 
mass function and the metric functions) also indicated a 
second-order convergence rate. 



V. NON-LINEAR COLLAPSE ON MINKOWSKI 

In this section we consider the case Mq = 0, Q = 0, 
namely, the collapse of the self-gravitating scalar field 
over a Minkowski spacetime, leading to the formation of 
a Schwarzschild-like black hole. Our initial data corre- 
spond to a compact sinusoidal ingoing pulse, as described 
in Section II. Here, we take v\ = 6, V% = 16, r$ = 6, and 
r„o = — 1/4 (corresponding to Mq = 0). The final mass of 
the black hole is then determined by the pulse amplitude 
A. In what follows we present the results of a numeri- 
cal simulation with A — 0.4, leading to a final black-hole 
mass Mf = 3.54. (Hereafter, we denote by Mf the final 
mass of the black hole.) 




10° 10 1 10 2 10 3 

N 

FIG. 6. The ! p -norms of the constraints C u and C v along an 
outgoing null ray, as functions of the grid parameter N. The 
cases a,b, and c refer to the h, h, and Zoo-norms, respectively, 
for C u , and cases d,e, and f refer to the h, h, and Zoo-norms, 
respectively, for C v . The numerical data are represented by 
circles, and the straight lines between the circles are linear 
interpolations of the data. 

Figure [?] displays the Bondi mass of the created black 
hole as a function of the retarded time w e .0 The Bondi 
mass decreases with u e , due to the escape of scattered 
energy to null infinity. The late-time decrease of the 
mass corresponding to the power-law tail of the scattered 
scalar field is too small to be observed in this figure (the 
numerical drift shown in Fig. || is also unobservable in 
the scale used in Fig. 0). 

The nonlinearity of the spacetime dynamics is best rep- 



**From the slopes of the curves displayed in Fig. |6| we can 
estimate the convergence rate of the code to be around 1.9, 
with variations of typical order 0.1. We stress, however, that 
these numbers would depend on the method employed for 
evaluating the convergence rate. 



''Strictly speaking, u e is not well-define here, as the space- 
time is dynamical and differs from Schwarzschild. In the 
asymptotic region r ^> Mf, however, the geometry becomes 
asymptotically Minkowski, and we can define u e with respect 
to this asymptotic region. Namely, along a ray v = const in 
this range, u e is linear with r. 
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resented by the evolution of the mass function along the 
AH (which is just twice the value of r there). Figure || 
shows this mass as a function of v. The mass function 
grows rapidly, until it approaches a saturation value. In 
this case, too, the mass-increase at late time due to the 
power-law tail, and the numerical mass drift, are unseen. 
The final black-hole mass can be deduced from either 
the flat l&rge-v portion of the graph in Fig. || or the flat 
large-Me portion in Fig. |?] — these two numbers agree, as 
they should. 



3.556 
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FIG. 7. The Bondi mass as a function of retarded time 
u e (calibrated such that u e — on u — 0), for N = 40. 
The mass is displayed along the ingoing null ray v = 10 6 M/, 
representing null infinity. 
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v 

FIG. 8. Black-hole mass determined from the AH radius 
vs. v. At advanced times earlier than v « 14 the domain of 
integration, u < it/, does not intersect the AH. 



The stability and accuracy of our code is demonstrated 
in Fig. which displays the scalar-field's QN ringing 



along the horizon for N = 5, 10, 20, and 40: The four 
graphs are indistinguishable in this Figure. In addition, 
we also determined the QN ringing frequency, and com- 
pared it with the linear analysis value |t7|]. This com- 
parison is hard, as we have only a few oscillations before 
the power-law tails start to dominate. In addition, the 
numerical mass drift (see above) complicates the com- 
parison between the numerical results and the theoretical 
prediction. (However, the mass drift can be controlled by 
the grid parameter N. Note that the physical mass in- 
crease due to scalar-field absorption is negligible at late 
times.) From our numerical data we find that the QN 
frequency is a = 0.032 — 0.026 i. The real part of a was 
calculated from the two nodes in Fig. ^ corresponding to 
a full wavelength, and the imaginary part from the two lo- 
cal extrema between them. The theoretical value for the 
least damped mode with I = is <j t h = 0.031 — 0.029 i 
(recall that here Mf = 3.54). The sources for the de- 
viation are the (numerical) drift in the mass, the effect 
of the other I = modes and the power-law tails, and 
the inability to use values from many cycles. However, 
our numerically obtained value is remarkably close to the 
linear analysis value. 
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FIG. 9. Quasi-normal ringing at the horizon as a function 
of v e . Recall that Mf ~ 3.54, which explains the relatively 
large value of v e in which the ringing takes place. We used 
here four different values of N — 5,10,20, and 40 — but the four 
graphs are indistinguishable in this figure. 

Figure |l^ shows the late-time behavior of $ in the 
three asymptotic regions: (a) - at fixed r, with t 3> Mf 
(we take t — (u e + u e )/2), (b) - at future null infinity 
(represented here by Vf = 10 6 M/), for u e s> Mf, and 
(c) - at the horizon, with v 3> Mf. This figure clearly 
demonstrates both the QN ringing and the power-law 
tails, in all three asymptotic regions. 

The determination of the asymptotic behavior at null 
infinity poses a special difficulty: We cannot integrate up 
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to v — oo proper. (An attempt to compactify the coordi- 
nate v will not solve this problem, as it would lead to a 
divergence of / at null infinity.) We therefore represent 
null infinity by a large (buy yet finite) value, v = Vf. 
This "null-infinity approximation" is only valid as long 
as Vf 3> u e . Thus, the determination of the late-time 
behavior at null infinity clearly demands huge values of 
Vf, in order to satisfy Vf 3> u e Mf. In the simula- 
tions described in this paper, we used Vf = 10 6 M/ to 
represent null infinity. In order to enable the integration 
to such large v values within a reasonable computation 
time, we used the following procedure: Let us denote 
by u e f the maximal value of u e in the desired presenta- 
tion of the late-time null-infinity behavior (in Fig. [l^, 
u e f = 10 4 ). After integrating up to a value of v which 
corresponds to v e = u e f, we chop the vectors u(I) and 
hi(I) at a value of u which corresponds to u e — u e f. The 
last point is now located at r > 2Mf. When we continue 
the integration to larger v values, the minimal value of 
r, r m i n (v) — r(u e — u e f,v), increases very rapidly and 
approaches large values (of order v). We can therefore in- 
crease dv accordingly, say, dv(v) ~ r m i n (v) /(ION). This 
allows us to integrate up to e.g., v = 10 10 within a very 
short integration time. (Practically, we change dv in dis- 
crete values of v, e.g. once in each decade.) 
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FIG. 10. The late-time behavior of the scalar field in the 
three limits (a),(b), and (c). Case (a): This graph displays 
$ at r = 2.3M/, as a function of t. Case (b): $ along null 
infinity (represented by v = 10 6 M/) versus retarded time u e , 
calibrated such that u e = on u = 0. Case (c): <& along 
the EH, as a function of v e . (The amplitude in this case 
was divided by 100, so that it will not overlap with the other 
graphs). In all three cases, the QN ringing and the power-law 
tails are seen clearly. We used here N = 20. 

As was mentioned above, one of our goals is to evalu- 
ate the power-law indices in the case of nonlinear collapse 
and to compare them to the predictions of the linear the- 
ory. Since our initial data correspond to I = and to 



zero initial static moment, the linear perturbation anal- 
ysis would predict (negative) power indices 3, 2, and 3 
in the cases (a), (b) and (c), correspondingly. In general, 
the slopes of the straight sections in the three graphs 
shown in Fig. |l^ appear to agree with these predicted 
values. However, the standard best-fit method is not 
so useful in this case for a precise determination of the 
numerically-computed indices, due to the following rea- 
son. Consider, for example, the late-time behavior at 
the horizon. According to the linear theory, it should 
be dominated by v~ 3 . However, this dominant term is 
"contaminated" by higher-order terms in 1/v, whose ef- 
fect become larger as v decreases Jl8| . Assume now that 
we use the standard best-fit method (applied to a finite 
interval v' < v < Vf) to determine the deviation of the 
power-law index from its predicted value. As it turns 
out, the computed deviation will be dominated in this 
case by the "higher-order contamination." This contam- 
ination effect, in turn, will depend in an arbitrary way 
on the choice of the parameter v' . In order to remove 
this arbitrariness, we introduce the notion of local power 
index, define by — v<$>, v /<f>. (For the other asymptotic 
regions, v is to be replaced by t or u e , accordingly.) 
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FIG. 11. Local index of the 'tails' for case (a): along 
r = const = 2.3M/, as a function of t. The local power is 
2.98 ± 0.01. We used here N = 20. 

The local power index for the three asymptotic regions 
is shown in Figs. [II], [jj], and [ll|. The agreement with the 
predictions of linear theory is remarkable. In principle, 
deviations from the precise integer index may result from 
three sources of errors: (i) the limited accuracy of the 
numerical simulation; (ii) the finiteness of the late-time 
domain covered by the numerics, i.e. the finiteness of t, u, 
and v in Figs. [II], [l2] and[l~5|, correspondingly (due to the 
"higher-orders contamination" , the precise integer index 
is expected only at infinitely-late time); (iii) in case (b) 
(i.e. at null infinity), the finiteness of the final value v = 



10 



v f taken to represent null infinity is also a possible source 
of error. In the numerical simulations presented here, we 
find that the deviation is related primarily to source (ii): 
We used a sufficiently large N, and a sufficiently large vt 
in case (b), so sources (i) and (iii) are insignificant. 



of the local power index. 

For comparison, we quote here the values obtained in 
previous nonlinear numerical analyses for the power-law 
indices. In case (a): 2.63-2.74 f| and 3.38 @; m case 
(c): 3.06 0; No nonlinear results where obtained so far 
for case (b). 
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FIG. 12. Local index of the 'tails' for case (b): along 
v = 10 6 M/ (representing future null infinity), as a function of 
u e . The local power is 2.002 ± 0.003. We used here N = 20. 
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FIG. 13. Local index of the 'tails' for case (c): along the 
AH, as a function of v e . The local power is 2.99 ± 0.02. We 
used here N = 20. 

Our results for the local index (at maximal i, u, or v) 
are: 

Case (a): 2.98 ± 0.01 (instead of 3), 
Case (b): 2.002 ± 0.003 (instead of 2), 
Case (c): 2.99 ± 0.02 (instead of 3). 
The error bar represents the numerical "noise" , produced 
primarily by the numerical differentiation of $ with re- 
spect to u, u, or t, which is inherent to the computation 



VI. NON-LINEAR COLLAPSE ON A CHARGED 
BACKGROUND 

In order to study the nonlinear dynamics of charged 
black holes, we consider here the gravitational collapse 
of the self-gravitating (neutral) scalar field over a pre- 
existing charged background (a RN geometry). The 
model and initial-value set-up are as explained in Sec- 
tion II. We take here an initial mass Mo = 1, and a 
charge Q — 0.95. (We found similar results for other val- 
ues of Q < 1.) We now take v\ — 6,t>2 = 16, r = 6, 
and r U Q w —0.1729. As before, we take a scalar-field 
amplitude A = 0.4. The black-hole mass then increases 
to Mf = 3.87 during the collapse. Figure |lj shows the 
value of r at the AH vs. v. The rapid increase of the 
horizon's area indicates strong nonlinear spacetime dy- 
namics. The two-stage increase of r reflects the structure 
of the scalar-field pulse: The latter has a maximum at 
about v = 11, and the vanishing of $ jt , implies M„ = 
there. (A similar behavior is observed in Fig. [?].) 
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FIG. 14. Value of r at the AH as a function of v. (At early 
values of v our numerical domain of integration u t < u < Uf 
does not intersect the AH.) We used here N = 20. 

According to the predictions of the linearized the- 
ory E3, the late-time behavior in the three asymptotic 
regions (a,b,c) in a (non-extreme) charged black hole 
should be similar to the uncharged case - namely, QN 
ringing followed by inverse power-law tails with the same 
indices as in the uncharged case. Figure Il5| displays the 
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late-time behavior of the scalar held for the three asymp- 
totic regions (a), (b), and (c). Again, both the QN ring- 
ing and the power-law tails are seen very clearly in all 
three asymptotic regions. 



Case (b): 1.996 ± 0.001 (instead of 2), 
Case (c): 2.99 ± 0.02 (instead of 3). 
These results are in excellent agreement with the predic- 
tions of the linear theory. 
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FIG. 15. Amplitude of the scalar field for the three cases 
(a), (b), and (c) for the nonlinear collapse on a charged RN 
background. Case (a): along r = const = 2.3Af/, as a func- 
tion of t. Case (b): along v = 10 6 Af/, representing future null 
infinity, as a function of u e (calibrated such that u e — on 
u — 0). Case (c): along the horizon, as a function of v e . The 
amplitude for case (c) is divided by 100 to avoid overlap of 
the graphs. We used here N = 20. 
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FIG. 16. Local power as a function of t for case (a): along 
constant value of r (r = 2.3M/). The value of the local power 
approaches 2.99 ± 0.01. We used here N = 20. 

Figures [l6|, [I?], and [l8] show the local power index for 
the three asymptotic regions. Our results for the local 
power index (at maximal t, u, or v) are: 
Case (a): 2.99 ± 0.01 (instead of 3), 




FIG. 17. Local power as a function of u e for case (b): along 
future null infinity, represented by v — 10 6 M/. The value 
of the local power approaches 1.996 ± 0.001. We used here 
N = 20. 
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FIG. 18. Local power as a function of v e for case (c): along 
the EH. The value of the local power approaches 2.99 ± 0.02. 
We used here N = 20. 



VII. CONCLUSIONS 

We developed a numerical scheme for the integration 
of the spherically-symmetric nonlinear Einstein-Maxwell- 
Klein-Gordon field equations. Our scheme is based on 
free evolution in double-null coordinates. This scheme is 
stable and accurate, it is capable of running to arbitrarily 
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late times, and it can handle black holes while avoiding 
singularities. 

We used this numerical code to study the gravitational 
collapse of a spherically-symmetric, self-gravitating, 
minimally-coupled scalar field to form a black hole (or the 
collapse of such a scalar field over a pre-existing charged 
background). Our numerical simulations demonstrate 
both the quasi- normal modes and the power-law tails, 
in all three late-time asymptotic regions: at a constant r 
(with large t) , along future null infinity (at large u) , and 
along the event horizon (at large v). The accuracy of 
our numerical scheme, its ability to run forever, and the 
method of calculating local power indices, allowed us to 
evaluate the power-law indices with an accuracy better 
than all previous estimates. 

Our results confirm that the predictions of the linear 
theory for the late-time behavior of perturbations out- 
side the black hole hold also for fully nonlinear collapse. 
(This observation is not surprising — in a sense, it is a 
manifestation of the principle that black holes have no 
hair.) In particular, in all three late-time asymptotic re- 
gions, the power-law indices approach asymptotically the 
integer values predicted by the linear perturbation anal- 
ysis. This agreement of the late-time nonlinear dynamics 
and the linear perturbation theory was already demon- 
strated by GPP H in the uncharged case (see also jOl ) . 
Here we demonstrate it for the charged case as well. 

The simulations presented here were restricted to the 
external part of the black hole and the neighborhood 
of the event horizon. One of our main motivations 
in this project, however, was to develop the numeri- 
cal tools which will allow the investigation of the inner 
part of black holes. We are currently using this numer- 
ical scheme to study the evolution of the geometry and 
the scalar field near the spacetime singularity inside a 
charged black hole. 

ACKNOWLEDGMENT 

This research was supported in part by the Israeli Sci- 
ence Foundation administered by the Israel Academy of 
Sciences and Humanities. 

APPENDIX A 

At issue here is the evolution of the entities at the left- 
hand sides of Eqs. (H|J), which represent the violation 
of the constraint equations. These entities are denoted 
in Ref . [^2| by E\ and E2 , respectively (in Section IV we 
denote these entities by C' u and C' v ). GP argue that, as 
a consequence of the precise evolution equations, E\ and 
E2 will grow exponentially with t along lines r = const. 
According to GP, this exponential divergence represents 
an inevitable numerical instability of the free-evolution 



scheme. We do not accept this conclusion, and claim 
that (i) under the precise evolution equations E\ and E2 
do not grow exponentially; rather, they are essentially 
conserved (in a sense which will be explained below); (ii) 
in our numerical free-evolution scheme, E\ (but not E2) 
will grow exponentially, but this growth is a consequence 
of the numerical error in the integration of the evolution 
equations, not of the equations themselves, (iii) This di- 
vergence of E\ does not indicate a numerical instability; 
Rather, it reflects the passive exponential growth of typi- 
cal gauge-dependent entities like r u along lines r = const 
in the Schwarzschild geometry. 

To verify point (i), assume that the evolution equa- 
tions are precisely satisfied. Then, Eq. (7) in Ref. JlJ] 
reads = —(r tV /r)E\. E2 will satisfy an analogous 
equation: E2, u — ~(r, u /r)E2- It then follows that the 
entity rE\ is conserved along lines u — const, and simi- 
larly rE2 is conserved along lines v — const. Therefore, 
an exponential growth of E\ or E2 along lines r = const 
is ruled outpj The exponential divergence of the linear 
metric perturbations (denoted £ and n in Ref. |l^]) found 
by GP must therefore be a gauge mode. In fact, the in- 
finitesimal coordinate transformation u — > u + du (e.g., 
with fixed infinitesimal du) yields a nonvanishing £ (in 
Ref. £ represents the linear perturbation in r 2 ), given 
by 

£ = — (r 2 ), u du — —2r r yU du . 

For any Kruskal-like u, at a fixed r, r u grows like 
exp(t/4Af) at f > M. It then follows that along any 
line r = const, at large t, £ will exhibit this exp(i/4M) 
divergence. This fits very well with the rate of diver- 
gence, exp(0.24i/M), found numerically by GP. But this 
is, of course, a gauge mode, which does not indicate a 
violation of the Einstein equations. 

Note also that the analytic derivation of the exponen- 
tial growth in Ref. Jl^] is based on the "mode ansatz" 
approximation. This approximation may only be valid 
if the mode's wave-number k is sufficiently large. The 
diverging modes found by GP do not satisfy this con- 
dition, however. This may explain the discrepancy of 
the value 0.32 M _1 predicted by GP compared with the 
above theoretical value 1/(4M) (which is also confirmed 
numerically by GP to a good accuracy). 

In our numerical tests, we found that along lines 
r = const, E2 is roughly preserved, while E\ grows like 
exp(i/2M). From the above discussion it is obvious that 
this exponential growth of E\ must be a consequence 



" Such an exponential divergence would demand that, on 
the initial hypersurface, E\ diverges exponentially with u e , 
and E2 diverges (almost) exponentially with v e , which is an 
unreasonable situation. 
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of the violation of the evolution equations, due to nu- 
merical errors. Later we shall give a more explicit ex- 
planation for this behavior. The crucial point is, how- 
ever, that the exponential growth of E\ does not indi- 
cate an exponential growth of the intrinsic local error: 
The entity E\ (like E-i) is not an appropriate error in- 
dicator, because it is not a gauge-invariant entity. [In a 
transformation u — > it' (it), E% is multiplied by the factor 
(du/du 1 ) 2 ] In order to extract from E\ the informa- 
tion about the intrinsic local error, we must construct 
a gauge-invariant entity from it. A convenient choice is 
the gauge-invariant entity C u = i?i/(r„) 2 . This entity 
indeed remains roughly constant along lines r = const 
and along the EH (see Fig. |J). The behavior of the 
other error indicators (e.g. the mass parameter in Figs, 
[l] and ||) also indicate stability: None of the invariant 
entities exhibit an exponential growth of error. 

We still need to explain why the numerical errors in 
the integration of the evolution equations results in an 
exponential growth of E\ . If there were no numerical er- 
rors in the integration, then, along a line r = const = 
r', Ei would approach (at large t) a constant value, 
Ei(uh,Vi) r(uh,Vi)/r' . Correspondingly, C u would decay 
like l/(r.„) 2 oc exp(— t/2M). However, the numerical er- 
ror in the integration of the evolution equations provide 
a (roughly) constant source of error in the evolution of 
Cj^. The combination of the dynamical tendency to ex- 
ponential decay and the constant source of error results 
in a finite saturation value (proportional to the numer- 
ical error)f**|. In turn, this implies that E\ grows like 
exp(i/2A/). 

Finally, we emphasize again that despite of the above 
discussion, whenever the domain of integration includes 
the EH, a naive attempt to use a free-evolution scheme in 
double-null coordinates, without a grid refinement, will 
inevitably result in a numerical instability due to the ex- 
ponential divergence of r.„ (see Section III). This might 
be the reason for the failure of previous attempts to use 
the free-evolution scheme. The grid refinement (in our 
case, the point-splitting procedure) is a necessary ingre- 
dient in such a numerical scheme. 
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^ To understand why this source of error is roughly inde- 
pendent of t, recall that (i) in our integration scheme, due to 
the dynamical grid refinement, the intrinsic error production 
rate is essentially independent of t (or u), and (ii) C u is an 
authentic indicator of the intrinsic local error. 

*** Phenomenologically, we may represent the situation by 
the simple differential equation dC u / dt = — C U /(2M) +S(r), 
where S(r) is the numerical source term. C u then approaches 
the asymptotic value 2MS(r). 
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